library(ggplot2)

Assume we have yield (or other agronomic assessment) data in spatial format, such as yield monitor data. For example, from SDSU

raw.dat <- read.csv("./data/sdsu/Swest_Corn2013-Clean.csv", header = TRUE)
summary(raw.dat)
##        ID             OID          Longitude        Latitude    
##  Min.   :    1   Min.   :    1   Min.   :-96.9   Min.   :43.23  
##  1st Qu.: 2905   1st Qu.: 8715   1st Qu.:-96.9   1st Qu.:43.23  
##  Median : 7054   Median :17428   Median :-96.9   Median :43.23  
##  Mean   : 8303   Mean   :17428   Mean   :-96.9   Mean   :43.23  
##  3rd Qu.:12993   3rd Qu.:26142   3rd Qu.:-96.9   3rd Qu.:43.23  
##  Max.   :21707   Max.   :34856   Max.   :-96.9   Max.   :43.24  
##            Field          Dataset          Product          Swath      
##  West_Field_SW:34856   Min.   :-113892   0533HR:34856   Min.   : 4.99  
##                        1st Qu.:-113892                  1st Qu.:14.99  
##                        Median :-113892                  Median :14.99  
##                        Mean   :-113892                  Mean   :14.92  
##                        3rd Qu.:-113892                  3rd Qu.:14.99  
##                        Max.   :-113892                  Max.   :14.99  
##     Distance        Duration     Track         Elevation    Areacount 
##  Min.   :0.820   Min.   :1   Min.   :  0.0   Min.   :1241   On:34856  
##  1st Qu.:4.530   1st Qu.:1   1st Qu.:  0.0   1st Qu.:1252             
##  Median :4.950   Median :1   Median :148.3   Median :1255             
##  Mean   :4.875   Mean   :1   Mean   :100.6   Mean   :1255             
##  3rd Qu.:5.250   3rd Qu.:1   3rd Qu.:180.0   3rd Qu.:1258             
##  Max.   :7.190   Max.   :1   Max.   :359.5   Max.   :1266             
##         Time          Yoffset     PassNum         CropFlw     
##  10/2/2013:34856   Min.   :0   Min.   : 1.00   Min.   : 0.40  
##                    1st Qu.:0   1st Qu.:12.00   1st Qu.:18.10  
##                    Median :0   Median :31.00   Median :20.60  
##                    Mean   :0   Mean   :34.24   Mean   :20.34  
##                    3rd Qu.:0   3rd Qu.:55.00   3rd Qu.:23.00  
##                    Max.   :0   Max.   :73.00   Max.   :35.60  
##     Moisture       YldMassDry        YldMassWet        YldVolDry     
##  Min.   : 0.00   Min.   :  299.9   Min.   :  307.1   Min.   :  5.35  
##  1st Qu.:16.87   1st Qu.:10866.3   1st Qu.:11173.1   1st Qu.:194.04  
##  Median :17.25   Median :11974.6   Median :12305.4   Median :213.83  
##  Mean   :17.28   Mean   :11837.6   Mean   :12165.9   Mean   :211.39  
##  3rd Qu.:17.63   3rd Qu.:13059.2   3rd Qu.:13424.0   3rd Qu.:233.20  
##  Max.   :24.00   Max.   :22241.7   Max.   :22736.6   Max.   :397.17  
##    YldVolWet          Speed            Prod          CropFlwV      
##  Min.   :  5.48   Min.   :0.560   Min.   :1.020   Min.   :  25.71  
##  1st Qu.:199.52   1st Qu.:3.090   1st Qu.:5.570   1st Qu.:1163.57  
##  Median :219.74   Median :3.380   Median :6.140   Median :1324.28  
##  Mean   :217.25   Mean   :3.324   Mean   :6.005   Mean   :1307.62  
##  3rd Qu.:239.71   3rd Qu.:3.580   3rd Qu.:6.460   3rd Qu.:1478.57  
##  Max.   :406.01   Max.   :4.900   Max.   :8.900   Max.   :2288.57  
##         Date      
##  10/2/2013:34856  
##                   
##                   
##                   
##                   
##

In this case, we’ll work with dry mass

hist(raw.dat$YldMassDry,main="YldVolDry")

ggplot(raw.dat, aes(Longitude, Latitude)) + geom_point(aes(colour = PassNum),size = 1)

We discard endrows for simplicity.

northBorder <- 43.23575
southBorder <- 43.22975
eastBorder <- -96.897
westBorder <- -96.900
min(raw.dat$Longitude)
## [1] -96.90052
max(raw.dat$Longitude)
## [1] -96.89669
trimmed.dat <- subset(raw.dat,raw.dat$Latitude>=southBorder)
trimmed.dat <- subset(trimmed.dat,trimmed.dat$Latitude<=northBorder)
trimmed.dat <- subset(trimmed.dat,trimmed.dat$Longitude<=eastBorder)
trimmed.dat <- subset(trimmed.dat,trimmed.dat$Longitude>=westBorder)

ggplot(trimmed.dat, aes(Longitude, Latitude)) + geom_point(aes(colour = YldVolDry),size = 1)

Now, let’s map a trial into one corner of the field. We will find this easier if we convert to meters. This will require an approximation, see http://stackoverflow.com/questions/639695/how-to-convert-latitude-or-longitude-to-meters

trimmed.dat$LonM <- trimmed.dat$Longitude - min(trimmed.dat$Longitude)
trimmed.dat$LatM <- trimmed.dat$Latitude - min(trimmed.dat$Latitude)
latMid <- (min(trimmed.dat$Latitude) + max(trimmed.dat$Latitude))/2
m_per_deg_lat = 111132.954 - 559.822 * cos( 2.0 * latMid ) + 1.175 * cos( 4.0 * latMid)
m_per_deg_lon = (3.14159265359/180 ) * 6367449 * cos ( latMid )
trimmed.dat$LonM <- trimmed.dat$LonM*m_per_deg_lon
trimmed.dat$LatM <- trimmed.dat$LatM*m_per_deg_lat
ggplot(trimmed.dat, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 1)

max(trimmed.dat$LonM)
## [1] 243.7949
max(trimmed.dat$LatM)
## [1] 666.539

So, let’s have a trial with plot dimensions of 2x4 m with 1 m buffer in between

plotDim <- c(2,4)
rowBuffer <- 1

We’ll start with Roy Scott’s 1993 plan

Scott.plan <- data.frame(
  row=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3),
  col=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
  trt=c(1, 5, 7, 2, 6, 8, 3, 10, 4, 9, 11, 12, 3, 14, 13, 1, 16, 2, 15, 4, 1, 19, 1, 2, 21, 18, 3, 20, 4, 17)
)
trial.dimensions <- function(plan,plot.dim,row.buffer) {
  rows <- max(plan$row)-min(plan$row)+1
  cols <-  max(plan$col)-min(plan$col)+1
  plot.width <- plot.dim[1]*cols + (cols-1)*row.buffer
  plot.height <- plot.dim[2]*rows + (rows-1)*row.buffer
  return(c(plot.width,plot.height))
}
trial.dim <- trial.dimensions(Scott.plan,plotDim,rowBuffer)

Add 3 meters to each side and subset our yield data

trial.dat <- subset(trimmed.dat,trimmed.dat$LonM<(trial.dim[1]+6))
trial.dat <- subset(trial.dat,trial.dat$LatM<(trial.dim[2]+6))
ggplot(trial.dat, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

Find points at the center of plots

Scott.plan$LonM <- 0
Scott.plan$LatM <- 0
start.point <- c(3,3)
half.width <- plotDim[1]/2
half.heigth <- plotDim[2]/2

for(idx in 1:dim(Scott.plan)[1]) {
  row <- Scott.plan$row[idx]
  col <- Scott.plan$col[idx]
  Scott.plan$LonM[idx] <- start.point[1] + half.width + (col-1)*plotDim[1] + (col-1)*rowBuffer
  Scott.plan$LatM[idx] <- start.point[2] + half.heigth + (row-1)*plotDim[2] + (row-1)*rowBuffer
}
ggplot(Scott.plan, aes(LonM, LatM)) + geom_point(aes(colour = row),size = 1)

Next, we krig the trial data

library(gstat)
sample.var <- variogram(YldVolDry~1, 
                        locations=~LonM+LatM, 
                        data=trimmed.dat)
plot(sample.var)

sample.vgm <- fit.variogram(sample.var, vgm("Exp"))
sample.vgm
##   model    psill   range
## 1   Nug 577.3469  0.0000
## 2   Exp 192.0298 23.6442
sample.krig <- krige(id="YldVolDry", 
                     formula=YldVolDry~1, 
                     data = trial.dat, 
                     newdata = Scott.plan, 
                     model = sample.vgm,
  locations=~LonM+LatM)
## [using ordinary kriging]
sample.krig
##    LonM LatM YldVolDry.pred YldVolDry.var
## 1     4    5       219.7577      624.6749
## 2     7    5       220.4189      622.1267
## 3    10    5       222.1915      622.5545
## 4    13    5       223.7153      624.6057
## 5    16    5       225.1562      622.4185
## 6    19    5       225.2944      622.7788
## 7    22    5       223.2662      622.0513
## 8    25    5       218.5649      621.1611
## 9    28    5       213.9409      622.8205
## 10   31    5       208.6830      623.9615
## 11    4   10       220.6846      622.4553
## 12    7   10       222.9458      620.2207
## 13   10   10       222.8527      620.7529
## 14   13   10       222.6355      622.1328
## 15   16   10       222.2083      620.1043
## 16   19   10       220.3901      621.0249
## 17   22   10       218.8548      619.9384
## 18   25   10       217.7635      619.5662
## 19   28   10       214.1294      620.8849
## 20   31   10       209.2803      621.3399
## 21    4   15       225.6756      624.3709
## 22    7   15       226.7191      622.4411
## 23   10   15       225.5107      622.8796
## 24   13   15       222.8781      624.1540
## 25   16   15       219.9744      621.7598
## 26   19   15       218.3965      622.3589
## 27   22   15       216.9719      622.1420
## 28   25   15       215.4486      621.5204
## 29   28   15       212.9803      622.8373
## 30   31   15       208.9109      623.9846
Scott.plan$YldVolDry <- sample.krig$YldVolDry.pred
ggplot(Scott.plan, aes(LonM, LatM)) + geom_point(aes(colour = row),size = 4)

pooled.dat <- data.frame(
  LonM = c(trial.dat$LonM,Scott.plan$LonM),
  LatM = c(trial.dat$LatM,Scott.plan$LatM),
  YldVolDry = c(trial.dat$YldVolDry,Scott.plan$YldVolDry)  
)
ggplot(pooled.dat, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

library(lme4)
## Loading required package: Matrix
summary(aov(YldVolDry ~ trt + row, data=Scott.plan))
##             Df Sum Sq Mean Sq F value Pr(>F)
## trt          1   19.0  19.044   0.703  0.409
## row          1    0.1   0.086   0.003  0.955
## Residuals   27  731.2  27.080
gy.lmer <- lmer(YldVolDry ~ trt + (1 | row), data=Scott.plan)
anova(gy.lmer)
## Analysis of Variance Table
##     Df Sum Sq Mean Sq F value
## trt  1 19.044  19.044  0.7292
print(gy.lmer, corr=FALSE) 
## Linear mixed model fit by REML ['lmerMod']
## Formula: YldVolDry ~ trt + (1 | row)
##    Data: Scott.plan
## REML criterion at convergence: 181.3479
## Random effects:
##  Groups   Name        Std.Dev.
##  row      (Intercept) 0.00    
##  Residual             5.11    
## Number of obs: 30, groups:  row, 3
## Fixed Effects:
## (Intercept)          trt  
##    220.5747      -0.1232
aug16.plan <- data.frame(
row=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8),
col=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)
)

Create a function to overlay a plan on a yield map

superimpose.plan <- function(plan,map.data,start.point,plot.dim=c(1,1),row.buffer=0,sample.vgm=NULL) {
  
  trial.dim <- trial.dimensions(plan,plotDim,row.buffer)
  
  trial.dat <- subset(map.data,map.data$LonM<(start.point[1]+trial.dim[1]+3))
  trial.dat <- subset(trial.dat,trial.dat$LatM<(start.point[2]+trial.dim[2]+3))

  trial.dat <- subset(trial.dat,trial.dat$LonM>=(start.point[1]-3))
  trial.dat <- subset(trial.dat,trial.dat$LatM>=(start.point[2]-3))
  
  plan$LonM <- 0
  plan$LatM <- 0
  
  half.width <- plot.dim[1]/2
  half.heigth <- plot.dim[2]/2
  
  for(idx in 1:dim(plan)[1]) {
    row <- plan$row[idx]
    col <- plan$col[idx]
    plan$LonM[idx] <- start.point[1] + half.width + (col-1)*plot.dim[1] + (col-1)*row.buffer
    plan$LatM[idx] <- start.point[2] + half.heigth + (row-1)*plot.dim[2] + (row-1)*row.buffer
  }
  if(is.null(sample.vgm)) {
    sample.var <- variogram(YldVolDry~1, 
                        locations=~LonM+LatM, 
                        data=map.data)
    sample.vgm <- fit.variogram(sample.var, vgm("Exp"))
  }
    
  sample.krig <- krige(id="YldVolDry", 
                     formula=YldVolDry~1, 
                     data = trial.dat, 
                     newdata = plan, 
                     model = sample.vgm,
  locations=~LonM+LatM)
  
  plan$YldVolDry <- sample.krig$YldVolDry.pred
  
  return(list(
    trial.dim=trial.dim,
    plan=plan,
    trial=trial.dat,
    krig=sample.krig,
    pooled = data.frame(
      LonM = c(trial.dat$LonM,plan$LonM),
      LatM = c(trial.dat$LatM,plan$LatM),
      YldVolDry = c(trial.dat$YldVolDry,plan$YldVolDry)  
    )
  ))
}
plan16sw <- superimpose.plan(aug16.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan16sw$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan16sw$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan16sw$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

plan16sw$trial.dim
## [1] 47 39
#plan16sw$plan
plan16se <- superimpose.plan(aug16.plan,
                           map.data=trimmed.dat,
                           start.point=c(190,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan16se$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan16se$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan16se$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

#plan16se$plan
plan16nw <- superimpose.plan(aug16.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,621),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan16nw$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan16nw$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan16nw$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

#plan16nw$plan
plan16ne <- superimpose.plan(aug16.plan,
                           map.data=trimmed.dat,
                           start.point=c(190,621),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan16ne$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan16ne$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan16ne$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

#plan16ne$plan
aug8.plan <- data.frame(
row=c(1,1,1,1,1,1,1,1,
      2,2,2,2,2,2,2,2,
      3,3,3,3,3,3,3,3,
      4,4,4,4,4,4,4,4,
      5,5,5,5,5,5,5,5,
      6,6,6,6,6,6,6,6,
      7,7,7,7,7,7,7,7,
      8,8,8,8,8,8,8,8,
      9,9,9,9,9,9,9,9,
      10,10,10,10,10,10,10,10,
      11,11,11,11,11,11,11,11,
      12,12,12,12,12,12,12,12,
      13,13,13,13,13,13,13,13,
      14,14,14,14,14,14,14,14,
      15,15,15,15,15,15,15,15,
      16,16,16,16,16,16,16,16),
col=c(1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1,
      1,2,3,4,5,6,7,8,
      8,7,6,5,4,3,2,1)
)
plan8sw <- superimpose.plan(aug8.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan8sw$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan8sw$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan8sw$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

plan8sw$trial.dim
## [1] 23 79
#plan8sw$plan
plan8se <- superimpose.plan(aug8.plan,
                           map.data=trimmed.dat,
                           start.point=c(218,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan8se$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan8se$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan8se$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

plan8se$trial.dim
## [1] 23 79
#plan8se$plan
plan8nw <- superimpose.plan(aug8.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,581),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan8nw$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan8nw$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan8nw$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

plan8nw$trial.dim
## [1] 23 79
#plan8nw$plan
plan8ne <- superimpose.plan(aug8.plan,
                           map.data=trimmed.dat,
                           start.point=c(218,581),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan8ne$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan8ne$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan8ne$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

plan8ne$trial.dim
## [1] 23 79
#plan8ne$plan
aug4.plan <- data.frame(
row=c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10,11,11,11,11,12,12,12,12,13,13,13,13,14,14,14,14,15,15,15,15,16,16,16,16,17,17,17,17,18,18,18,18,19,19,19,19,20,20,20,20,21,21,21,21,22,22,22,22,23,23,23,23,24,24,24,24,25,25,25,25,26,26,26,26,27,27,27,27,28,28,28,28,29,29,29,29,30,30,30,30,31,31,31,31,32,32,32,32),
col=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4)
)
plan4 <- superimpose.plan(aug4.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan4$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan4$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan4$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

#plan4$plan

Spatially balanced plans

plan8x14.plan <- data.frame(
row=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8,8),
col=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14),
vanEs14x8=c(1,12,14,6,5,8,2,7,10,11,13,4,9,3,5,7,1,9,8,4,11,13,6,12,10,2,3,14,7,10,11,14,9,1,6,12,13,3,4,5,2,8,6,4,10,1,3,5,8,14,11,13,9,12,7,2,11,8,3,12,1,7,4,6,2,9,14,10,13,5,13,1,2,11,4,14,7,3,5,6,12,9,8,10,12,11,5,4,10,2,1,9,14,8,3,7,13,6,4,14,8,7,13,12,10,1,3,5,2,11,6,9)
)

plan8x14.plan$blk <- plan8x14.plan$row
plan8x14.plan$trt <- plan8x14.plan$vanEs14x8

plan8x14sw <- superimpose.plan(plan8x14.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
ggplot(plan8x14sw$trial, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 6)

ggplot(plan8x14sw$plan, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

ggplot(plan8x14sw$pooled, aes(LonM, LatM)) + geom_point(aes(colour = YldVolDry),size = 4)

aov.tbl <- summary(aov(YldVolDry ~ as.factor(trt)+as.factor(blk),data=plan8x14sw$plan))
aov.tbl
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## as.factor(trt) 13    574    44.1   0.899    0.557    
## as.factor(blk)  7   3663   523.3  10.663 9.81e-10 ***
## Residuals      91   4466    49.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
trial.dim <- trial.dimensions(plan8x14.plan,plot.dim=c(2,4),row.buffer=1)

trial.width <- trial.dim[1]
trial.height <- trial.dim[2]


TrtMS <- c()
TrtP <- c()
RepMS <- c()
RepP <- c()
ResMS <- c()
ResP <- c()
Row <- c()
Col <- c()
        
rightBorder <- max(trimmed.dat$LonM) - (trial.width+6)
topBorder <- max(trimmed.dat$LatM) - (trial.height+6)

row=1
col=1
atColEnd=FALSE
atRowEnd=FALSE
while(!atColEnd) {
  currentRow <- 3+(row-1)*trial.height + (row-1)*6
  if(currentRow < topBorder) {
    while(!atRowEnd) {
      currentCol <- 3+(col-1)*trial.width + (col-1)*6
      if(currentCol<rightBorder) {
        corner <- c(currentCol,currentRow)
print(corner)
        #overlay and analyze
        currentPlan <- superimpose.plan(plan8x14.plan,
                           map.data=trimmed.dat,
                           start.point=corner,
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
        aov.tbl <- summary(aov(YldVolDry ~ as.factor(vanEs14x8)+as.factor(row),data=currentPlan$plan))
        
        TrtMS <- c(TrtMS,aov.tbl[[1]][1,3])
        TrtP <- c(TrtP,aov.tbl[[1]][1,5])
        RepMS <- c(RepMS,aov.tbl[[1]][2,3])
        RepP <- c(RepP,aov.tbl[[1]][2,5])
        ResMS <- c(ResMS,aov.tbl[[1]][3,3])
        ResP <- c(ResP,aov.tbl[[1]][3,5])
        Row <- c(Row,corner[1])
        Col <- c(Col,corner[2])
        
        col=col+1  
      } else {
        atRowEnd=TRUE
      }
    }
    col=1
    row=row+1
    atRowEnd=FALSE
  } else {
    atColEnd =TRUE
  }
}
## [1] 3 3
## [using ordinary kriging]
## [1] 50  3
## [using ordinary kriging]
## [1] 97  3
## [using ordinary kriging]
## [1] 144   3
## [using ordinary kriging]
## [1] 191   3
## [using ordinary kriging]
## [1]  3 48
## [using ordinary kriging]
## [1] 50 48
## [using ordinary kriging]
## [1] 97 48
## [using ordinary kriging]
## [1] 144  48
## [using ordinary kriging]
## [1] 191  48
## [using ordinary kriging]
## [1]  3 93
## [using ordinary kriging]
## [1] 50 93
## [using ordinary kriging]
## [1] 97 93
## [using ordinary kriging]
## [1] 144  93
## [using ordinary kriging]
## [1] 191  93
## [using ordinary kriging]
## [1]   3 138
## [using ordinary kriging]
## [1]  50 138
## [using ordinary kriging]
## [1]  97 138
## [using ordinary kriging]
## [1] 144 138
## [using ordinary kriging]
## [1] 191 138
## [using ordinary kriging]
## [1]   3 183
## [using ordinary kriging]
## [1]  50 183
## [using ordinary kriging]
## [1]  97 183
## [using ordinary kriging]
## [1] 144 183
## [using ordinary kriging]
## [1] 191 183
## [using ordinary kriging]
## [1]   3 228
## [using ordinary kriging]
## [1]  50 228
## [using ordinary kriging]
## [1]  97 228
## [using ordinary kriging]
## [1] 144 228
## [using ordinary kriging]
## [1] 191 228
## [using ordinary kriging]
## [1]   3 273
## [using ordinary kriging]
## [1]  50 273
## [using ordinary kriging]
## [1]  97 273
## [using ordinary kriging]
## [1] 144 273
## [using ordinary kriging]
## [1] 191 273
## [using ordinary kriging]
## [1]   3 318
## [using ordinary kriging]
## [1]  50 318
## [using ordinary kriging]
## [1]  97 318
## [using ordinary kriging]
## [1] 144 318
## [using ordinary kriging]
## [1] 191 318
## [using ordinary kriging]
## [1]   3 363
## [using ordinary kriging]
## [1]  50 363
## [using ordinary kriging]
## [1]  97 363
## [using ordinary kriging]
## [1] 144 363
## [using ordinary kriging]
## [1] 191 363
## [using ordinary kriging]
## [1]   3 408
## [using ordinary kriging]
## [1]  50 408
## [using ordinary kriging]
## [1]  97 408
## [using ordinary kriging]
## [1] 144 408
## [using ordinary kriging]
## [1] 191 408
## [using ordinary kriging]
## [1]   3 453
## [using ordinary kriging]
## [1]  50 453
## [using ordinary kriging]
## [1]  97 453
## [using ordinary kriging]
## [1] 144 453
## [using ordinary kriging]
## [1] 191 453
## [using ordinary kriging]
## [1]   3 498
## [using ordinary kriging]
## [1]  50 498
## [using ordinary kriging]
## [1]  97 498
## [using ordinary kriging]
## [1] 144 498
## [using ordinary kriging]
## [1] 191 498
## [using ordinary kriging]
## [1]   3 543
## [using ordinary kriging]
## [1]  50 543
## [using ordinary kriging]
## [1]  97 543
## [using ordinary kriging]
## [1] 144 543
## [using ordinary kriging]
## [1] 191 543
## [using ordinary kriging]
## [1]   3 588
## [using ordinary kriging]
## [1]  50 588
## [using ordinary kriging]
## [1]  97 588
## [using ordinary kriging]
## [1] 144 588
## [using ordinary kriging]
## [1] 191 588
## [using ordinary kriging]
hist(TrtMS)

hist(TrtP)

hist(RepMS)

hist(RepP)

hist(ResMS)

sim.dat <- data.frame(
          TrtMS = TrtMS,
        TrtP = TrtP,
        RepMS = RepMS,
        RepP = RepP,
        ResMS = ResMS,
        Row = Row,
        Col = Col
)

ggplot(sim.dat, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(sim.dat, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(sim.dat, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

And a function to simulate a plan over a field

simulate.plan <- function(plan,field,plot.dim=c(1,1),row.buffer=0,sample.vgm=NULL) {
  trial.dim <- trial.dimensions(plan,plot.dim=plot.dim,row.buffer=row.buffer)

  trial.width <- trial.dim[1]
  trial.height <- trial.dim[2]

  if(is.null(sample.vgm)) {
    sample.var <- variogram(YldVolDry~1, 
                        locations=~LonM+LatM, 
                        data=map.data)
    sample.vgm <- fit.variogram(sample.var, vgm("Exp"))
  }
    
  TrtMS <- c()
  TrtDF <- c()
  TrtP <- c()
  RepMS <- c()
  RepDF <- c()
  RepP <- c()
  ResMS <- c()
  ResP <- c()
  ResDF <- c()
  Row <- c()
  Col <- c()
        
  rightBorder <- max(field$LonM) - (trial.width+6)
  topBorder <- max(field$LatM) - (trial.height+6)

  row=1
  col=1
  atColEnd=FALSE
  atRowEnd=FALSE
  while(!atColEnd) {
    currentRow <- 3+(row-1)*trial.height + (row-1)*6
    if(currentRow < topBorder) {
      while(!atRowEnd) {
        currentCol <- 3+(col-1)*trial.width + (col-1)*6
        if(currentCol<rightBorder) {
          corner <- c(currentCol,currentRow)
          #overlay and analyze
          currentPlan <- superimpose.plan(plan,
                           map.data=field,
                           start.point=corner,
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
          aov.tbl <- summary(aov(YldVolDry ~ as.factor(trt)+as.factor(blk),data=currentPlan$plan))
          
          
          TrtDF <- c(TrtDF,aov.tbl[[1]][1,1])
          TrtMS <- c(TrtMS,aov.tbl[[1]][1,3])
          TrtP <- c(TrtP,aov.tbl[[1]][1,5])
          RepDF <- c(RepDF,aov.tbl[[1]][2,1])
          RepMS <- c(RepMS,aov.tbl[[1]][2,3])
          RepP <- c(RepP,aov.tbl[[1]][2,5])
          ResDF <- c(ResDF,aov.tbl[[1]][3,1])
          ResMS <- c(ResMS,aov.tbl[[1]][3,3])
          ResP <- c(ResP,aov.tbl[[1]][3,5])
          Row <- c(Row,corner[1])
          Col <- c(Col,corner[2])
        
          col=col+1  
        } else {
          atRowEnd=TRUE
        }
      }
      col=1
      row=row+1
      atRowEnd=FALSE
    } else {
      atColEnd =TRUE
    }
  }
  sim.dat <- data.frame(
               TrtDF = TrtDF,
               TrtMS = TrtMS,
               TrtP = TrtP,
               RepDF = RepDF,
               RepMS = RepMS,
               RepP = RepP,
               ResDF = ResDF,
               ResMS = ResMS,
               Row = Row,
               Col = Col
            )
  return(sim.dat)
}

Block Size 14

Stagger 0

aug14.plan <- data.frame(
row=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8,8),
col=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14,1,2,3,4,5,6,7,8,9,10,11,12,13,14),
trt=c(81,1,2,3,82,4,5,6,83,7,8,9,84,10,84,11,12,13,81,14,15,16,83,17,18,19,82,20,83,21,22,23,81,24,25,26,82,27,28,29,84,30,82,31,32,33,81,34,35,36,84,37,38,39,83,40,83,41,42,43,84,44,45,46,81,47,48,49,82,50,84,51,52,53,83,54,55,56,82,57,58,59,81,60,81,61,62,63,83,64,65,66,84,67,68,69,82,70,81,71,72,73,82,74,75,76,84,77,78,79,83,80)
)
aug14.plan$blk <- aug14.plan$row

trial.dimensions(aug14.plan,plot.dim=c(2,4),row.buffer=1)
## [1] 41 39
stag0 <- simulate.plan(aug14.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag0, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag0, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag0, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

Stagger 1

aug14.plan$trt <- c(81,1,2,3,82,4,5,6,83,7,8,9,84,10,20,83,11,12,13,84,14,15,16,81,17,18,19,82,83,30,84,21,22,23,81,24,25,26,82,27,28,29,39,81,40,83,31,32,33,82,34,35,36,84,37,38,48,49,82,50,84,41,42,43,83,44,45,46,81,47,57,58,59,84,60,83,51,52,53,81,54,55,56,82,83,67,68,69,81,70,84,61,62,63,82,64,65,66,76,81,77,78,79,84,80,82,71,72,73,83,74,75)
stag1 <- simulate.plan(aug14.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag1, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag1, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag1, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

Stagger 2

aug14.plan$trt <- c(81,1,2,3,82,4,5,6,83,7,8,9,84,10,83,20,84,11,12,13,81,14,15,16,82,17,18,19,28,29,84,30,81,21,22,23,82,24,25,26,83,27,84,37,38,39,81,40,82,31,32,33,83,34,35,36,45,46,81,47,48,49,82,50,83,41,42,43,84,44,83,54,55,56,82,57,58,59,81,60,84,51,52,53,62,63,81,64,65,66,84,67,68,69,83,70,82,61,83,71,72,73,84,74,75,76,82,77,78,79,81,80)
stag2 <- simulate.plan(aug14.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag2, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag2, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag2, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

aug14_7.plan <- data.frame(
row=c(1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,3,3,4,4,4,4,4,4,4,5,5,5,5,5,5,5,6,6,6,6,6,6,6,7,7,7,7,7,7,7,8,8,8,8,8,8,8,9,9,9,9,9,9,9,10,10,10,10,10,10,10,11,11,11,11,11,11,11,12,12,12,12,12,12,12,13,13,13,13,13,13,13,14,14,14,14,14,14,14,15,15,15,15,15,15,15,16,16,16,16,16,16,16),
col=c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6,7),
trt=c(81,1,2,3,82,4,5,83,6,7,8,84,9,10,81,11,12,13,84,14,15,82,16,17,18,83,19,20,82,21,22,23,84,24,25,83,26,27,28,81,29,30,83,31,32,33,81,34,35,84,36,37,38,82,39,40,81,41,42,43,82,44,45,84,46,47,48,83,49,50,82,51,52,53,84,54,55,83,56,57,58,81,59,60,82,61,62,63,81,64,65,84,66,67,68,83,69,70,83,71,72,73,82,74,75,84,76,77,78,81,79,80)
)

aug14_7.plan$blk <- ceiling(aug14_7.plan$row/2)

stag0b <- simulate.plan(aug14_7.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag0b, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag0b, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag0b, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

aug14_7.plan$trt <- c(81,1,2,3,82,4,5,10,83,6,7,8,84,9,14,15,82,11,12,13,83,81,19,20,84,16,17,18,23,82,24,25,83,21,22,27,28,84,29,30,81,26,31,32,33,81,34,35,82,83,36,37,38,84,39,40,45,82,41,42,43,81,44,49,50,84,46,47,48,83,82,54,55,83,51,52,53,58,84,59,60,81,56,57,62,63,81,64,65,82,61,66,67,68,84,69,70,83,84,71,72,73,83,74,75,80,82,76,77,78,81,79)

stag1b <- simulate.plan(aug14_7.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag1b, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag1b, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag1b, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

aug14_7.plan$trt <- c(81,1,2,3,82,4,5,9,10,83,6,7,8,84,13,81,14,15,82,11,12,16,17,18,83,19,20,84,25,82,21,22,23,84,24,83,29,30,81,26,27,28,32,33,81,34,35,82,31,84,36,37,38,83,39,40,44,45,83,41,42,43,81,48,84,49,50,82,46,47,51,52,53,82,54,55,81,60,83,56,57,58,84,59,83,64,65,84,61,62,63,67,68,81,69,70,82,66,84,71,72,73,82,74,75,79,80,81,76,77,78,83)

stag2b <- simulate.plan(aug14_7.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag2b, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag2b, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag2b, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

stag0$stagger <- "14.0"
stag1$stagger <- "14.1"
stag2$stagger <- "14.2"
stag0b$stagger <- "7.0"
stag1b$stagger <- "7.1"
stag2b$stagger <- "7.2"

stag.dat <- rbind(stag0,stag1,stag2,stag0b,stag1b,stag2b)

stag.dat$stagger <- as.factor(stag.dat$stagger)
summary(aov(TrtMS ~ stagger,data=stag.dat))
##              Df  Sum Sq Mean Sq F value Pr(>F)
## stagger       5   13448    2690   0.574   0.72
## Residuals   393 1840492    4683
plot(TrtMS ~ stagger,data=stag.dat)

summary(aov(TrtP ~ stagger,data=stag.dat))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## stagger       5  14.10  2.8205   43.62 <2e-16 ***
## Residuals   393  25.41  0.0647                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(TrtP ~ stagger,data=stag.dat)

summary(aov(RepMS ~ stagger,data=stag.dat))
##              Df   Sum Sq Mean Sq F value   Pr(>F)    
## stagger       5  2964297  592859   10.66 1.27e-09 ***
## Residuals   393 21853553   55607                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(RepMS ~ stagger,data=stag.dat)

summary(aov(RepP ~ stagger,data=stag.dat))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## stagger       5  16.56   3.312   39.48 <2e-16 ***
## Residuals   393  32.97   0.084                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(RepP ~ stagger,data=stag.dat)

summary(aov(ResMS ~ stagger,data=stag.dat))
##              Df  Sum Sq Mean Sq F value Pr(>F)    
## stagger       5  128074   25615    8.25  2e-07 ***
## Residuals   393 1220227    3105                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ResMS ~ stagger,data=stag.dat)

plot(TrtDF ~ stagger,data=stag.dat)

plot(RepDF ~ stagger,data=stag.dat)

plan16x7.plan <- plan8x14.plan
plan16x7.plan$row=aug14_7.plan$row
plan16x7.plan$col=aug14_7.plan$col


plan16x7.plan$blk <- ceiling(plan16x7.plan$row/2)

stag.vanes7 <- simulate.plan(plan16x7.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag.vanes7, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag.vanes7, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag.vanes7, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

hist(stag.vanes7$TrtMS)

hist(stag.vanes7$TrtP)

hist(stag.vanes7$RepMS)

hist(stag.vanes7$RepP)

hist(stag.vanes7$ResMS)

max(stag0$Row)
## [1] 191
max(stag0$Col)
## [1] 588
max(stag0b$Row)
## [1] 211
max(stag0b$Col)
## [1] 513

8 rows by 14 columns vanEs14x8 aug80v140n, etc.

max(stag.vanes7$Row)
## [1] 211
max(stag.vanes7$Col)
## [1] 513
plan8x14sw <- superimpose.plan(plan8x14.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan8x14sw$plan

plan8x14se <- superimpose.plan(plan8x14.plan,
                           map.data=trimmed.dat,
                           start.point=c(211,513),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan8x14se$plan

plan8x14nw <- superimpose.plan(plan8x14.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan8x14nw$plan

plan8x14ne <- superimpose.plan(plan8x14.plan,
                           map.data=trimmed.dat,
                           start.point=c(211,513),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan8x14ne$plan
stag.vanes7 <- simulate.plan(plan16x7.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
max(stag.vanes7$Row)
## [1] 211
max(stag.vanes7$Col)
## [1] 513

16 rows by 7 columns

aug80v70n, etc.

plan16x7sw <- superimpose.plan(plan16x7.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan16x7sw$plan

plan16x7nw <- superimpose.plan(plan16x7.plan,
                           map.data=trimmed.dat,
                           start.point=c(211,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan16x7nw$plan

plan16x7se <- superimpose.plan(plan16x7.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,513),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan16x7se$plan

plan16x7ne <- superimpose.plan(plan16x7.plan,
                           map.data=trimmed.dat,
                           start.point=c(211,513),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan16x7ne$plan
plan32x4.plan <- data.frame(
row=c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,9,9,9,9,10,10,10,10,11,11,11,11,12,12,13,13,13,13,14,14,14,14,15,15,15,15,16,16,17,17,17,17,18,18,18,18,19,19,19,19,20,20,21,21,21,21,22,22,22,22,23,23,23,23,24,24,25,25,25,25,26,26,26,26,27,27,27,27,28,28,29,29,29,29,30,30,30,30,31,31,31,31,32,32),
col=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2,1,2,3,4,1,2,3,4,1,2,3,4,1,2),
trt=c(1,12,14,6,5,8,2,7,10,11,13,4,9,3,5,7,1,9,8,4,11,13,6,12,10,2,3,14,7,10,11,14,9,1,6,12,13,3,4,5,2,8,6,4,10,1,3,5,8,14,11,13,9,12,7,2,11,8,3,12,1,7,4,6,2,9,14,10,13,5,13,1,2,11,4,14,7,3,5,6,12,9,8,10,12,11,5,4,10,2,1,9,14,8,3,7,13,6,4,14,8,7,13,12,10,1,3,5,2,11,9,6))
plan32x4.plan$blk <- ceiling(plan32x4.plan$row/4)
stag.vanes4 <- simulate.plan(plan32x4.plan,trimmed.dat,plot.dim=c(2,4),row.buffer=1,sample.vgm=sample.vgm)
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
ggplot(stag.vanes4, aes(Row, Col)) + geom_point(aes(colour = RepP),size = 8)

ggplot(stag.vanes4, aes(Row, Col)) + geom_point(aes(colour = ResMS),size = 8)

ggplot(stag.vanes4, aes(Row, Col)) + geom_point(aes(colour = TrtP),size = 8)

hist(stag.vanes4$TrtMS)

hist(stag.vanes4$TrtP)

hist(stag.vanes4$RepMS)

hist(stag.vanes4$RepP)

hist(stag.vanes4$ResMS)

max(stag.vanes4$Row)
## [1] 224
max(stag.vanes4$Col)
## [1] 498
plan32x4sw <- superimpose.plan(plan32x4.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan32x4sw$plan



plan32x4nw <- superimpose.plan(plan32x4.plan,
                           map.data=trimmed.dat,
                           start.point=c(3,498),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan32x4nw$plan

plan32x4ne <- superimpose.plan(plan32x4.plan,
                           map.data=trimmed.dat,
                           start.point=c(224,498),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan32x4ne$plan

plan32x4se <- superimpose.plan(plan32x4.plan,
                           map.data=trimmed.dat,
                           start.point=c(224,3),
                           plot.dim=c(2,4),
                           row.buffer=1,sample.vgm=sample.vgm
                           )
## [using ordinary kriging]
#plan32x4se$plan